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ABSTRACT 


The behavior of Fe-Cr alloys under irradiation is in part controlled by the 
characteristics of point defects generated by high energy collision. Radiation enhanced 
diffusion and radiation induced precipitation are among the mechanisms that lead to 
changes in the microstructure under irradiation, and are thus controlling effects such as 
swelling and a’ precipitation. Point defects in Fe-Cr alloys are diverse in nature due to 
their interaction with a variety of local solute configurations. Ab initio results indicate 
that the magnetic structure of the alloy is critical in determining its energetics. The ability 
to model these properties with classic potentials is still to be proven. In this work a 
detailed comparison between ab initio and classic values of a variety of point defects 
configurations is performed, testing in this way the extent to which classic potentials can 
be reliably used for radiation damage studies, and evaluating the dependence of point 


defect formation energies on Cr concentration. 
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I. INTRODUCTION 


The Fe-Cr alloys are considered for use as structural materials for Generation IV 
advanced nuclear designs. High doses of neutron irradiation cause vacancies and 
interstitials to form in the steel matrix. The most accurate approach for evaluating 
formation energies of point defects is numerical integration based on the first principles 
of quantum mechanics. However, due to the computational requirements posed by this 
method, it is prohibitive to use it for large scale simulations. It is therefore necessary to 
provide an approximation of ab initio results using Molecular Dynamics and Monte Carlo 
methodologies coupled with an adequately descriptive semi-empirical many body 
potential for alloys. The potential evaluated in this work is based on combining recent ab 
initio results with thermodynamic properties to describe formation energy as a function of 
local composition. This thesis seeks to answer the following questions: 

1. Determine the extent to which this potential adequately approximates 
results obtained by ab initio calculations in pure elements. 

2. Determine how formation energies of point defects vary as a function 
of Cr concentration; determine their mutual relationships and 
stability in all possible configurations. 

In order to answer these questions, simulations were conducted using the 
Lawrence Livermore National Laboratory’s massively parallel super-computing 
hardware and a suite of publicly available and custom designed materials science 
software solutions. 

The calculations of this thesis show that reasonable predictions can be obtained 
using the potential described in this work. Formation energies of point defects in pure 
elements and as a function of Cr concentration and heretofore unexpected behavior of 


formation energies versus Cr concentration have been discovered. 
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Il. BACKGROUND 


A. HISTORY OF NUCLEAR POWER 


The use of nuclear power for energy generation traces its origins to the volatile 
times preceding World War II. The theoretical basis for fission was laid by Otto Hahn 
and Fritz Strassman in the German periodical Naturwissenschaften in January 1939, 
where they reported that an isotope of barium was produced by neutron bombardment of 
uranium. A knowledge network formed connecting Otto Robert Frisch and Lise Meitner, 
German scientists escaping Hitler’s rule to Denmark, with Niels Bohr who was on his 
way to visit Albert Einstein at Princeton. Once the information got to Princeton which 
had become a Mecca for theoretical physics when Dr. Einstein joined its Institute for 
Advanced Studies, it set off a “chain reaction” in the scientific community. By the end of 
the year over a hundred papers on the subject of fission had been published in Physical 
Review and the concept of a sustainable chain reaction was grounded in theory [1]. 

The engineering required to sustain and control a chain reaction was, at the time, 
far from trivial. Nonetheless, the first artificial nuclear reactor, Chicago Pile 1 (CP-1), 
was built and reached criticality on December 2", 1942 at the underground racquetball 


courts at the University of Chicago [2]. 






Figure 1. Drawing of Chicago Pile 1 [3] 


At the time there was no possibility of enriching the nuclear fuel, so natural 
uranium had to be used, which contains only about 0.7 % of the highly fissionable 


3 


uranium-235 isotope, with the remainder consisting of the more stable uranium-238. It 
was therefore necessary to use large quantities of highly purified graphite as neutron 
moderation material. The moderator was necessary to slow down the high energy 
neutrons produced by fission to an energy range where they would be much more easily 
absorbed by U-235 around .025 eV, roughly the same temperature as that of the 
surrounding material; hence the name “thermal.” [4] The fission cross section of U-235 


as a function of neutron energy is shown in Figure 2. 


NEUTRON CROSS-SECTIONS FOR FISSION OF URANIUM AND PLUTONIUM 
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Figure 2. Neutron Cross-Sections for Fission of Uranium and Plutonium [5] 


As we can see from Figure 2, neutrons produced by fission will fall within the 
energy range of .1 MeV to 5 MeV, but uranium-235 has the highest fission cross section 
at energies well below that value. The “moderator” is therefore used to slow down the 
fission-produced neutrons and increase their probability of causing fissions. The 
“moderator” is usually a highly purified substance, generally deuterium (in the form of 
heavy water), beryllium, or graphite for natural uranium reactors. These materials do not 
absorb neutrons easily (they have a low absorption cross section), but their nuclei are 
light enough that neutrons collide inelastically, transferring their kinetic energy to the 
moderator and therefore slowing down in the energy spectrum towards the energies 
which make them more likely to cause further fissions [4]. A diagram showing a typical 


moderated chain reaction is shown in Figure 3: 
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Figure 3. A Moderated Chain Reaction 


Once it became possible to enrich uranium, therefore increasing the percentage of 
uranium 235 in the fuel, it became possible to utilize light water for both moderation and 
cooling of large scale engineering systems. This advance generally heralded the arrival of 
the next generation of nuclear reactors. This generation of reactors was, for the most part, 
typified by pressurized water reactors (PWR) and boiling water reactors (BWR) in the 
West and the VVER/RBMK reactors in the former Soviet Union and satellite states. 
Canada, due to the abundance of natural uranium ore, opted for a heavy water reactor 
design (CANDU) which used natural rather than enriched uranium. All of these systems 
are considered to be second generation nuclear reactors, two common threads of which 
are the unique design and features of each reactor unit and their thermal-neutron based 


operating cycle [6]. 


B. FUTURE DIRECTIONS IN NUCLEAR POWER 


The newly emerging, third generation of reactor designs are equipped with 
advanced features such as safety systems incorporating passive energy dissipation or 
natural processes, simplifying their design and allowing them to cope with malfunctions 
without the need for operator action. Even though these designs are still confined to the 
thermal end of the fission cross section, third generation designs show some marked 
improvements when compared to their second generation cousins. Some of those are 


listed below: 


e astandardized design for each type to expedite licensing, reduce capital cost and 
reduce construction time, 

e asimpler and more rugged design, making them easier to operate and less 
vulnerable to operational upsets, 

e higher availability and longer operating life - typically 60 years, 

e reduced possibility of core melt accidents, 

e minimal effect on the environment, 


e higher burn-up to reduce fuel use and the amount of waste [7]. 


After the Three Mile Island accident, United States imposed a self-imposed 
moratorium on building new nuclear reactors. While third generation reactors have been 


built in Japan, Taiwan and Europe, this technology stagnated in the United States. 


In 2000, however, eleven world nations formed The Generation IV International 
Forum (GIF), a research and development consortium tasked with leading the way 
toward innovative nuclear energy systems. Based on eight far-ranging technology goals, 
Generation IV nuclear energy systems are aimed at achieving nuclear energy’s potential 


worldwide. The objective is a new generation of nuclear energy systems that: 


e advance nuclear safety; 
e address nuclear nonproliferation and physical protection issues; 
e are competitively priced, and 


e minimize waste and optimize natural resource utilization. 


Among the six most promising designs for Generation IV, three are thermal 
neutron spectrum systems (the Very-High-Temperature Reactor (VHTR) and 
Supercritical-Water-Cooled Reactor (SCWR)) with coolants and temperatures that enable 
hydrogen or electricity production with high efficiency, and three are fast neutron 
spectrum systems (the Gas-Cooled (GFR), the Lead-Cooled (LFR), and the Sodium- 
Cooled (SFR) fast reactors) that will enable more effective management of actinides 
through recycling of most components in the discharged fuel. [8] 

A significant change in Generation IV relative to previous reactor generations is 
the planned commercialization of fast spectrum systems. A fast spectrum neutron reactor 
is a reactor in which the chain reaction is sustained by fast neutrons without significant 
thermalization. Such a reactor needs no neutron moderator, but must generally use fuel 
that is relatively rich in fissile material when compared to that required for a thermal 
reactor. A diagram outlining the differences between fast and thermal spectrum processes 


is shown in Figure 4. 
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Figure 4. Overview of Differences Between Fast and Thermal Absorption Processes 


Because parasitic absorption in the moderator can result in a major loss of 
neutrons in a thermal reactor, a fast reactor has an inherently superior neutron economy; 
that is, there are excess neutrons not required to sustain the chain reaction. These 
neutrons can be used to produce extra fuel, as in the fast breeder reactor, or to transmute 
long half-life waste to less troublesome isotopes, such as in the Super Phénix reactor near 
Cadarache in France, or in some combination of these two purposes [9]. An overview of 


nuclear energy development is shown in Figure 5: 
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Figure 5. Overview of Reactor Technology Development [10] 


C. FAST SPECTRUM NUCLEAR REACTOR DESIGN FEATURES 


An example of a fourth generation system currently under development is a Lead- 
Cooled Fast Reactor (LFR), which features a fast-spectrum lead or lead/bismuth eutectic 
liquid metal-cooled reactor with a closed fuel cycle. The Pb coolant is a poor absorber of 
fast neutrons and this enables the realization of improved sustainability and fuel cycle 
benefits. Pb does not interact vigorously with air, water/steam, or carbon dioxide, 
eliminating concerns about exothermic reactions. It has a high boiling temperature, so the 
prospect of boiling or flashing of the ambient pressure coolant is realistically eliminated. 
The LFR is cooled by natural convection with a reactor outlet coolant temperature of 550 
°C, possibly ranging up to 800 °C with advanced materials. The higher temperature 


enables the production of hydrogen by thermochemical processes [11] 





Lead-Cooled Fast Reactor 








Figure 6. Lead-Cooled Fast Reactor Diagram [10] 


Additionally, due to its neutron economy, an LFR can be designed in such a 
manner as to behave like a converter or even a breeder reactor, therefore enabling the 
conversion of the uranium 238 isotope into fissionable plutonium 239. A neutron 
captured by uranium 238 results in the formation of uranium 239, which has a half life of 
about 23 minutes and decays into neptunium 239 through beta decay. Neptunium 239 has 
a half life of 2.4 days and then decays into plutonium 239, also through beta decay [12]. 
The added benefit of the high neutron flux and fast spectrum core design is the ability of 
the LFR to serve in the actinide management role. [13] 

Since the fission of fuel materials initially present in the core creates more new 
fuel (plutonium 239) from uranium 238, as a side product of fission, a properly designed 
reactor core can significantly extend the time between refueling, or indeed eliminate it 
altogether. An advanced reactor design that maximizes this property is the Small, Sealed, 
Transportable, Autonomous Reactor (SSTAR), currently under development by 


Lawrence Livermore and Argonne National Laboratories [11]. 
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D. SSTAR PROGRAM OVERVIEW 


The objective of the SSTAR program is to create a sealed reactor that can be 
delivered to a site, left to generate power for up to 30 years, and retrieved when its fuel is 
spent. The potential for nuclear proliferation would be minimized by sealing the cartridge 
core inside a tamper proof cask. The reactor would be monitored and operated 
autonomously with a satellite link to a national authority or international authorities 
overseeing the reactors. The system would include passive safety features such as a fast 
spectrum core with a strong reactivity feedback that enables autonomous load-following 
and provide passive power shutdown in case of a loss-of-coolant accident, as well as 
natural circulation flow of the Pb coolant. This is a critical issue, as the reactor’s primary 
market is thought to be in undeveloped or underdeveloped regions of the world, where it 
would operate with almost full autonomy. This marketing strategy also provides the 
reasoning for the small size of the reactor, since conventional nuclear stations typically 
produce about a gigawatt of electricity, which would overwhelm the distribution grid in a 


developing country, therefore wasting much of the installed power [14]. 
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Figure 7. Diagram of a Potential SSTAR Design [14] 
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The nominal SSTAR design (20 MWe/ 45 MWt) is 18 meters tall and 3.2 m in 
diameter, which enables it to be easily moved by barge or rail. The producer would be 
responsible for delivering the sealed unit by ship and truck and installing it at the 
operating location. At the end of its regular operating life, the old reactor would be 


removed for recycling or disposal, and a replacement system would be deployed. [14] 





Figure 8. SSTAR Deployment [14] 


SSTAR would be coupled with an S-CO, gas turbine Brayton cycle power 
converter, which enables the reactor to operate with a higher efficiency when compared 
to the traditional Rankine saturated steam cycle. This modification is enabled by the 
higher core temperatures in the LFR design, attainable with the primary Pb coolant and 
highly enriched transuranic (TRU) nitride fuel in '"N in a compact core. The temperatures 
in consideration are 650°C peak cladding temperature and 561°C core outlet temperature 
for a 405°C core inlet temperature [13]. However, the advantages of this design face a 
critical limitation on the lifetime of this and other Generation IV reactor designs from the 


material corrosion which results from radiation damage to structural materials 
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Ul. MOTIVATION 


A. RADIATION DAMAGE 


Fission-emitted neutrons carry with them up to 5 MeV of energy. Continuous 
neutron bombardment of structural materials exposed to such high energy neutron 
radiation results in introduction of defects within the material’s crystalline lattice. The 
inelastic collisions between the impacting neutrons and the structural atoms cause energy 
to be transferred to the atoms. Since the structural materials are, for the most part, 
comprised of stable elements, the system response to the higher energy state of the atom 
is its displacement from its original site within the lattice. This disorder within the perfect 
crystal lattice is referred to as a Frenkel disorder [15]. It results in the formation of a 
Frenkel pair of defects, where there is a single vacancy and a single interstitial within the 


structure [15]. 


Vacant 
lattice 
sité 





Figure 9. Vacancy Site in a Crystalline Structure [16] 


Vacancies are sites which would normally be occupied by an atom but which as a 
result of the displacement are unoccupied. If a neighboring atom moves to occupy the 
vacant site, the vacancy moves in the opposite direction to the site which used to be 
occupied by the moving atom. The stability of the surrounding crystal structure 
guarantees that the neighboring atoms will not simply collapse around the vacancy. In 
some materials, neighboring atoms actually move away from a vacancy, because they can 


better form bonds with atoms in the other directions [15]. 
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Interstitials are atoms which occupy a site in the crystal structure at which there is 
usually not an atom. They are generally high energy configurations. There are two basic 
kinds of interstitials: Intrinsic and extrinsic interstitials. Intrinsic interstitials are 
interstitial atoms of the same kind as the atoms of the crystal ("self-interstitials"). They 
are practically nonexistent in elemental crystals (that are found in all metals) with the 
significant exception of Si, where intrinsic interstitials play an important role in diffusion 
and microdefect formation. Extrinsic interstitials are interstitial atoms of a foreign 
(extrinsic) type, e.g. C in Fe or O in Si (“mixed interstitials”). They may diffuse 
directly through the lattice (1.e., without the help of vacancies) and play an important role 


in many technically relevant materials [15] 


Interstitial 
atom 





Figure 10. Interstitial Atom in a Crystal Lattice [16] 


Complexes can form between different kinds of point defects. For example, if a 
vacancy encounters an impurity, the two may bind together if the impurity is too large for 
the lattice. Interstitials can form 'split interstitial’ or 'dumbbell' structures where two 
atoms effectively share an atomic site, resulting in neither atom actually occupying the 
site [15]. These interstitial pairs can be oriented in three possible directions: 


e <100>: Atoms displaced in only in one axis as shown in Figure 11; 
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Figure 11. <100> Interstitial [17] 


e <110>: Atoms displaced in two axes as shown in Figure 12; and 











Figure 12. <110> Interstitial [17] 


e <111>: Atoms displaced in all three axes as shown in Figure 3-8. 
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Figure 13. <111> Interstitial [17] 


The result of a single neutron collision is not a single defect; on the contrary the 


interaction of a single 1 MeV neutron with an atomic nucleus transfers up to about 100 


15 


keV to the primary knock-on atom (PKA) of iron. Some of this recoil energy is lost to 
interactions with the electron cloud, resulting in a somewhat lower kinetic energy that is 
dissipated in atomic collisions. The PKA kinetic energy is then transferred by numerous 
subsequent collisions and resultant displacements, producing further generations of 
recoiling atoms at lower energies. The process terminates when the kinetic energy of the 
nth-generation of recoils falls below that needed to cause additional displacements [18]. 
Closely spaced interstitials and vacancies quickly recombine and only about one-third of 
the initial displacements survive. Typically, this leaves a vacancy-rich cascade core, 
surrounded by a shell of interstitials. The majority of interstitials quickly cluster to form 
small, disc-shaped features that are identical to small dislocation loops [19]. Along with 
interstitials, these loops are very mobile. Diffusion of individual interstitials and loops 
within the cascade region causes additional recombination prior to their rapid long-range 
migration (unless they are strongly trapped by other defects or solutes). Although they are 
less mobile than interstitials, vacancies also eventually diffuse. Through a series of local 
jumps, the vacancies and solutes in the cascade quickly begin to evolve to lower energy 
configurations, forming small, three-dimensional clusters, while others leave the cascade 
region [19]. The small clusters are unstable and can dissolve by vacancy emission. 
However, the small clusters also rapidly diffuse and coalesce with each other, forming 
larger nanovoids, which persist for much longer times. Solute atoms bind to the vacancies 
and segregate to clusters. The vacancy emission rate is lower from vacancy-solute cluster 
complexes. Small solute clusters remain after all the vacancy clusters have finally 
dissolved. Expressing damage exposure, or neutron dose, in terms of displacements-per- 
atom (dpa) partially accounts for the effect of the neutron energy spectrum on the 


generation of cascade defects and the net residual defect production scales with dpa [20]. 


B. MACRO-LEVEL EFFECTS OF RADIATION DAMAGE 


This time evolution of microscopic defects results in the following emergent 
macroscopic material behaviors: 
e Volumetric Swelling: irradiated material will display volumetric 


expansion in all axes. This results from lattice parameter elongation and 
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void formation [21]. An example of radiation induced swelling of 


autensitic steel is shown in Figure 14 
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Figure 14. Radiation Induced Swelling [21] 
Irradiation Creep: a _ time-dependent, constant rate mechanical 
deformation of material occurring slowly at stresses below the ultimate 
tensile stress. Vacancies cause biased absorption of interstitials at 
dislocations leading to net dislocation climb [22]. An example of 


irradiation creep is shown in Figure 15. 





Figure 15: Irradiation Creep of Nuclear Fuel Pin Cladding [23] 
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e Radiation hardening and embrittlement. This is a second order 
phenomenon which is caused by the slipping over the dislocation-linked 
clusters of point defects. This leads to decreased plasticity of the irradiated 
material and ultimately greater proclivity towards rapid catastrophic 
fracture generation. These effects become even more pronounced as the 


ambient temperature increases [18]. 
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Figure 16: Diagram of Material Hardening and Embrittlement Evolution [18] 


C. MATERIAL ISSUES FOR GENERATION IV NUCLEAR ENERGY 
SYSTEMS 


The fact that nuclear materials are subject to radiation damage is well known from 
the many reactor years of operation accumulated in operation of Generations I and II 
systems. However, when we compare the operating conditions present in early 
generations of nuclear reactors to those projected for Generation IV and other advanced 
nuclear applications, we see that the operating environment adversity significantly 


increases. 
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Fission Fission Fusion NASA 
(Gen. I) (Gen. IV) (Demo) space react. 


Structural alloy <300°C 500-1000°C 550-1000°C ~1000°C 
maximum temperature 








Max dose for core ~1 dpa ~30-100 dpa ~150 dpa ~10 dpa 
internal structures 





Max transmutation ~0.l appm | ~3-10 appm ~1500 appm ~] appm 
helium concentration (~10000 appm 

for SiC) 

Coolants H,O He, H,O, Pb- | He, Pb-Li, Li | Li, Na, or 
Bi, Na He-Xe 


Structural Materials Zircaloy, | Ferritic steel, Ferritic/ Nb-1Zr, Ta 
stainless SS, martensitic alloys, Mo 
steel superalloys, | steel, V alloy, alloys, 
C- composite | SiC composite | superalloys 




















Table 1. Operating Conditions for Advanced Nuclear Applications [24] 


This increase in adversity, as seen from Table 1, make the materials previously 
used for nuclear applications, namely austenitic steels (such as 316-L), no longer suitable 
for advanced applications. Therefore, ferritic-martensitic steels, such as HT-9, T-91 and 
EP823 have to be considered. The chemical composition of the ferritic-martensitic steels 


is compared to an example of autensitic steel in Table 2 


PEE PEPE PPP EP 
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Table 2. Chemical Composition of Potential Generation IV Structural Materials [25] 


We can see that a common thread among the materials considered is that they are 
all variants of iron-chromium (Fe-Cr) alloy, but the following are the practical reasons 


why these materials were selected: 
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1. Unlike stainless steel which starts to swell almost immediately after being 
irradiated, Fe-Cr alloys exhibit a delay in swelling, wherein they do not 
exhibit that property until receiving 100 dpa of radiation damage. 

2. Fe-Cr alloys perform well at high temperatures and pressures required for 
operation of advanced nuclear reactors, and they have favorable rates of 
corrosion when exposed to the elements intended for use as coolants (Pb, 
Bi, Na). 

3. Their welding properties are well known and there is experience with high 
temperature helium embrittlement [26]. 

However, despite all these favorable properties, the behavior of defect clustering 
resulting in irradiation creep, hardening and embrittlement is poorly understood in Fe-Cr 
alloys. Considering that these behaviors are known to occur in materials used in earlier 
reactor designs, which operate in much more benign conditions, it becomes imperative to 
thoroughly understand them. Additionally, when we take into account the intended life- 
time of a Generation IV system (around 60 reactor years) and that some designs (such as 
SSTAR) are intended to operate autonomously for extended periods of time, it is clear 
that our ability to model the time evolution of Fe-Cr alloys is critical for safe operation of 


these systems. 


D. MULTI-SCALE APPROACH TO MATERIAL SCIENCE MODELING 


Since radiation induced defects occur at the atomistic level, the most accurate 
method of calculating their energetic properties is to solve the Schrédinger equation (or 
Dirac equation, if relativistic effects are important) in its many-body electron form. 
However, this is an exceptionally difficult and computationally intensive process, as the 
currently “standard” model for condensed matter physics, the Density Functional Theory 
(DFT) using Local Density Approximation (LDA), is currently limited to 100-1000 
atoms. The largest DFT simulation to date is 1080 B atoms (with 3840 electrons 
simulated) on Lawrence Livermore National Laboratory’s 2000 CPU Linux cluster [23]. 
While the DFT approach is generally successful in predicting structures and macroscopic 
properties, it under-predicts band gap energies, over-predicts lattice parameters, and 
predicts wrong ground states for some magnetic systems (e.g., Fe). Nonetheless, for the 
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purposes of material and reactor engineering, requiring time evolution of macro scale 
behaviors such as void swelling, hardening, embrittlement, creep, stress corrosion 
cracking, the first-principles calculation is not suitable [27]. Therefore, in order to arrive 
at a successful time-evolved simulation of a finite structural element, several layers of 
approximation are necessary. 

The multi-scale modeling approach to material science is only natural, because 
the introduction of defects, their evolution and mutual interaction and ultimately their 
effect on the mechanical properties of the material all occur on very broad timescales, as 


can be seen in Figure 17. 
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Figure 17. Temporal and Spatial Placement of Phenomena [23] 


The modeling of these phenomena in a global sense is a chain of events in which the 
results of each previous step serve as the source of inputs and basis of comparison for the 
subsequent steps. Fitting the model at the next higher level of approximation is done by 
comparing its results for properties best modeled by the method on a finer scale. By using 
a bottom-up approach, based on first principles and built upon scale reversible models 
(when possible) starting from the sub-nanometer scale (where the building blocks of 
matter are established, hence providing material unity and technology integration), 
complete characterization of properties in materials and processes at different 
dimensional and time scales can be achieved. In addition, a sequence of experiments is 


designed to follow the modeling effort and validates the simulation result. A visual 
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description of the overall multi-scale modeling effort, complete with supporting 


experiments is shown in Figure 18. 
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Figure 18. Global View of Multi-scale Material Modeling Effort [28] 
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IV. THEORY AND METHODOLOGY FOR COMPUTATION 


A. INTRODUCTION 


Modern materials science evolved from two basic elements--metallurgy and 
alchemy--and has over the period of the last couple of decades greatly expanded its field 
of study. Its natural complexity had, for a very long time, kept its progress in check as 
the difficult physical concepts could only be applied to idealized, pure forms. 
Nonetheless, Moore’s Law compounding of computer power has recently enabled 
numerical simulation of increasingly complex systems. While the most exact calculations 
are based on first principles of quantum mechanics, alternatively referred to as ab initio 
calculations, their sheer complexity has yet to be conquered by Moore’s Law. The 
complexity of these calculations is clear to anyone who has ever taken an introductory 
course in quantum mechanics, which culminates in solving the Schrédinger Equation for 
a hydrogen atom. However, since most useful materials contains significantly more 
complex elements, the difficulty is compounded as it becomes necessary to solve 


Schrédinger’s many-electron equation in the following form: 


N 2 N 

AW=(T+V4U|¥= > -*Vv? + Y>V (7) + > U(F;, r;)|W=EW 

: eg (1) 
Where H is the electronic molecular Hamiltonian, N is the number of electrons and U is 
the electron-electron interaction. The operators T and U are so-called universal operators 
as they are the same for any system, while V is system dependent or non-universal. In 
essence, the greatest complexity in ab initio methods arises from the many-body 
interaction term U. While there are many clever methods for reducing the computational 
power required to solve this equation for a complex system, those remain a discussion for 
another day. A powerful approach has been developed in Ref. [29], that uses classic 
interatomic potentials and incorporates ab initio and experimental results, in order to 


enable the study of crystalline defects and their long range interactions on a much larger 


scale. 
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B. THE MANY-BODY POTENTIAL 


The models generally used for this work are colloquially known as ‘‘many-body”’ 
potentials, and can be globally grouped in two categories: the embedded-atom models 
and the second moment approximation [30]. However, the issue facing material modelers 
for advanced nuclear applications is that the overwhelming body of work in this field is 
focused on either pure elements or intermetallic compounds; only a few authors address 
concentrated alloys. 

The methodology used for this work was developed by Dr. Alfredo Caro, et al. 
of the Lawrence Livermore National Laboratory, first published in the Physical Review 
Letters article of August 2005 [29]. The objective of this methodology is to address 
arbitrarily complex systems of concentrated alloys with complex heat of formation. What 
makes this work especially appropriate is that this methodology is applied to Fe-Cr alloys 
which are of interest for advanced nuclear applications. 

The “many-body” potentials are based on the summation of atom energies 
which provides the total energy of the system. The atom energies are themselves 
composed of two contributions; namely, embedding and pair potential terms. For 


heteroatomic systems, such as binary alloys composed of elements A and B, this reads: 


N 
Ee = >| Fa Sone, (rs ) + 5D Va0, (rs) | 


jJ#i ~ j#l 


(2) 
where a and £ stand for elements A and B sitting at sites i and j within the crystalline 
lattice, F’s are the embedding functions for either type of element, and V’s and p’s are the 
pair potentials and densities between a-f pairs. The functions p,g and Vg therefore 
describe the properties of the alloy. The variety in expressions of embedding energies, 
densities, and pair potentials encompasses a great similitude of models. 

The greatest issue in application of this model to the case of Fe-Cr is that unlike 
other binary alloys, its formation energy is not a symmetric function. Instead, in the case 
of Fe-Cr this function is highly nonsymmetrical and it even changes sign at low Cr 


composition [35], as seen in figure 19. 
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Figure 19. Ab Initio Calculations of Fe-Cr Mixing Enthalpy [35] 


Therefore, in contrast to the binary alloys with symmetric formation energy 
such as Au-Ni and Fe-Cu, the standard approach using a cross pair potential term is not 
good enough to reproduce the properties of the systems. By using the potentials already 
described in the literature, adjusting the alloy term in equation (2) and focusing on 
nonlinearities built upon the pair potential cross term alone, this methodology manages to 
create a model which successfully departs from ideality and matches the heat of 


formation behavior as described by Par Olsson et al. in [35] and as shown in Figure 20. 
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Figure 20. Formation Energy of the Alloy as Predicted by the Potential Used in This 
Work [29] 
C, THE FREE ENERGY EXPRESSION 


The free energy model uses an effective representation of two pure element 
potentials as reported by Mendelev for Fe [36] and Wallenius for Cr [37], with 


normalized densities, which for a = A, B reads 
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Pa ™ Pl Cnew 
re is Fo(@a) meet el (Qieq) Co 


Vaal?) = Ve,a(r) + 2F°"(O%,)p2(r), 
scl (3) 


where the superscript ° stands for original, 7° ,, for the density on a lattice site at 


eq 
equilibrium Dae P.,(1;")], and the prime ’ for derivative. These transformations do not 


alter the properties of the pure elements but have the advantage of minimizing the 
contribution of the embedding term to the formation energy of the alloy, thus enabling 
unrelated pure element potentials to be combined in this alloy description, as shown by 
the free energy of a random solution alloy with composition x at temperature T in 


equation 4: 


2(x, T) = Sree(x, T) + 2mix(x, T) + Ag(x, T), 
(4) 


where rer is the compositional weighted free energy of the pure components, gmix is the 
free energy contribution from the entropy of mixing for a random alloy, and Ag is the 
excess Gibbs energy of mixing, which, when expressed by a Redlich-Kister expansion 
[38] reads 


Ag(x,T) =x(1 - yy EU = 2x)F, 
9=0 
’ (5) 


where L, is the p-th order binary interaction parameter. This parameter is also the object 
of two major simplifications in this methodology, as it is generally a function of 
temperature. Therefore, in order to reduce the complexity of the description we neglect 
the excess vibrational entropy and assume that the formation energy does not depend on 


T. This leads to a simplified version of equation (5). 


Ag(x, T) = AH(x) = x(1— x) > L,(1 — 2x). 
p=0 (6) 


The factors for the Redlich-Kister expansion come from the Par Olsson ab initio 


calculations [35] and are given in the table below. 
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Lo Ly L2 L3 L4 








0.41566 0.0814134 -0.0101899 0.267659 -0.248269 

















Table 3.Values of the Redlich-Kister expansion coefficients for Eq. (6) 


D. THE CROSS PAIR POTENTIAL 


The functional form of the cross potential is based on the analytic mode of the 
alloy in which the species that sits on site i can be either A or B, but both are embedded 


in the same average environment, as discussed by Ackland and Vitek [31]: 


Etand = x sy VaalTi;) + er ~ VpplT ij) 


+ 2x4Xp > Vaplrij) + x4F 4(P) + xpF (A) 


(7) 
where the p is defined as 
p =X, palrij) + xpLpzslrij) ne 
Therefore, the contribution of the embedding terms to the mixing energy, AE™” is 
AE™ = x,[F4(p) — Falp = 1)] 
+ xglFp(p) — Fp(p = 1] . 


Using a Taylor expansion of F around p=l and using equation (3), it becomes clear that 


this contribution is quadratic in (p - 1), and therefore small for small variations in p-1 


as seen in 


AEM = x, Fi (p = 1)(p — 1)? + xg Fi(p = 1p — 1). dé 
The specific potentials used for this study reduce the embedding term contribution 
to the formation energy down to ~1 meV/atom at x = 0.5, making it negligible in 
comparison to the target value of ~100 meV for Fe-Cr alloys [35]. This enables us to 
write down the energy contribution from the pair potential terms as follows, 
AEPr = x(1 — x){2uap — (ug + vp)}- 


(11) 
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It is therefore assumed that: a) the nonlinearity is dependent only on the pair 
potential and b) Vaz is a function of both x and r, separable as a product of h(x)uap(r), 


with the form, 


Vaplx, r) = ACxY{Vaa(r) + Ven (r)]. = 


The h(x) is also a 4" order Reidlich-Kister polynomial, whose factors are obtained 
through global minimization of equation (7) using Mathemathica™ , and listed in Table 


4. 





Ho hy H2 H3 hy 








0.883644 -0.059302 0.644634 -1.342524 0.918932 

















Table 4. 4" Order Polynomial Coefficient for h(x), as Obtained by Minimization 


This methodological approach yields a formation energy curve indistinguishable 
from the target function obtained by the ab initio calculations of Par Olsson [35], as seen 
in Figure 19. Therefore, it can be concluded that it is reasonable to utilize this potential 


for modeling the formation energies of point defects in Fe-Cr alloys. 


E. THE MOLECULAR DYNAMICS CODE 


According to Ohno et al. [42], the method of classical molecular dynamics (MD) 
was first proposed by Alder and Wainwright in 1957. The essence of the approach comes 
from two assumptions: a) that there is a set of defined interatomic force functions or a 
potential, as we have shown in the beginning part of this chapter; and b) the numerical 
ensemble in which the number of particles N, the temperature T and the chemical 
potential are constant. If the above two statements are correct, the equation of motion 


for the atoms is the usual Newtonian equation 


2 N 
I ee YT (13) 


mM; = ng i 
dt? 





where m; is the mass of the i-th particle, r; its coordinate and F; the force acting on it. The 
code solves this equation numerically and performs energy minimization of the ensemble, 


therefore achieving thermal equilibrium after a sufficient number of iterations. This in 
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turn provides us with the energetics of a sample in which all the atoms have achieved 
their least energetic state and are in a “relaxed” state. 

One of the greatest limitations of the MD approach has been in coping with the 
effects of finite system size and presence of surfaces. In order to reduce these effects, the 
code uses periodic boundary conditions. With this approach, all the particles are placed 
inside a box, or a unit cell. If the particle goes outside the cell, it is brought back in from 
the opposite side of the cell. While this approach works marvelously with perfect crystal 
systems, introducing a disturbance such as a point defect creates problems because the 
code is inclined to see a “mirror image” of the defect across the periodic boundary 
conditions, therefore unintentionally creating a cascading effect of mutual influence. 

One technique for minimizing the effects of finiteness of the simulated system 
is to sum the forces exerted on the i-th particle from all other particles inside a certain 
“sphere of influence” defined by the cut-off radius R, centered on the particle. Of course 
this summation must be performed even if the particles are not in the same unit cell, but 
in the image cells [43]. Ideally, having an adequately large sample size in which the 
sphere of influence radius around the i-th atom does not encompass any atoms affected 


by the presence of its mirror image solves this problem. 
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A. 


V. PROCEDURE AND RESULTS 


RESEARCH PLAN 


The thesis research was conducted according to the following schedule: 


Task 1: The Study of Defects in Pure Elements. The objective of this section of 
the thesis is to prove that the Caro potential [29] adequately approximates results 
obtained by ab initio calculations. This will be done by calculating formation 
energies of point defects in pure element samples and comparing those results 
against values available in literature, obtained through ab initio calculations. 

Task 2: Evaluating the Vacancy Energy of Formation. This section is aimed at 
evaluating the behavior of vacancy formation energy in samples with increasing 
chromium concentration. This will serve to define the specific energetics of this 
point defect for a range of Cr concentration in relation to a linear interpolation 
between vacancy formation energies in pure elements. 

Task 3: Obtaining the Formation Energies of Interstitials. In this section the 
research will focus on defining the behavior of formation energies as a function of 
Cr concentration for a variety of interstitial configurations, as well as exploring 


their relative relationships. 


Put together, the scope of these results will serve to prove the applicability of this 


type of a semi-empirical potential for approximating the results of ab initio 


calculations and provide static and dynamic calculations on a scale currently 


unobtainable by ab initio simulation. 
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B. CHARACTERIZING THE PERFORMANCE OF THE CODE: THE 
STUDY OF DEFECTS IN PURE ELEMENTS 


In order to benchmark the performance of the newly developed potential [29] 
against other available Fe and Cr potentials and methodologies, it was necessary to 
compare its prediction of formation energies for the 12 point defects types of interest: 

e Single Vacancy in a pure Fe sample 

e Single Vacancy in a pure Cr sample 

e Self-interstitial (Fe-Fe) <110> in a pure Cr sample 

e Self-interstitial (Fe-Fe) <111> in a pure Fe sample 

e Self-interstitial (Fe-Fe) <110> in a pure Fe sample 

e Self-interstitial (Fe-Fe) <111> in a pure Cr sample 

e Mixed interstitial (Fe-Cr) <110> in a pure Fe sample 

e Mixed interstitial (Fe-Cr) <111> ina pure Fe sample 

e Self-interstitial (Cr-Cr) <110> in a pure Cr sample 

e Self-interstitial (Cr-Cr) <111> in a pure Cr sample 

e Mixed interstitial (Fe-Cr) <110> in a pure Cr sample 

e Mixed interstitial (Fe-Cr) <111> ina pure Cr sample 
The benchmarking could only be performed at the endpoints of the Cr concentration 
diagram, since no other physically sound models for the Fe-Cr alloys could be found in 
literature. 

The simulation samples were prepared by manually modifying the 1024 atom, 
perfect crystal samples of pure Fe and pure Cr. For the vacancy cases, this modification 
consisted of deleting the atom located at coordinates (0,0,0) in both samples. For the 
interstitial cases, these modifications consisted of adding an additional atom (atom 1025) 
to the end of the sample and shifting the first atom on the list (atom 1). The location of 
atom 1025 was determined by adding 0.5 A to the coordinates in the axes where the 
displacement was to occur, i.e., in x and y axes for <110> type interstitials and in all three 
axes for <111> type interstitials. The value of 0.5 A was subtracted from the respective 
coordinates of atom 1, therefore creating the relevant interstitial pair. The value of 0.5 A 
was determined by trial and error, as initial attempts to use a larger initial interstitial 


spacing (0.9 A) produced results in which the interstitial pair would fly apart during 
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energy minimization. For self-interstitials, an atom of the same type would be inserted 
into the sample, while the mixed interstitial samples received an atom of the opposite 
type. 

The samples were then run, using LLNL’s Multiprogramatic Capability Cluster — 
MCR. This is a large (11.2 TF) tightly coupled Linux cluster with 1,152 nodes, each with 
two 2.4-GHz Pentium 4 Xeon processors and 4 GB of memory. MCR runs the LLNL 
CHAOS software environment, which incorporates the Red Hat Linux operating system 
[44]. The Molecular Dynamics (MD) code, Moldy [45], was set up for constant 
temperature of 0.01 K, no pressure, no chemical potential, with free sample volume, and 
a run of 3000 steps with a time step of 10° seconds. This number of steps was necessary 
in order to ensure that the sample had ample time to relax around the induced defect and 
that the final energy values converged. At the end of the run, the final total enthalpy 
values were extracted for each output and compared against those obtained from running 
the perfect crystal samples under the same conditions. The table listing the final values of 


total enthalpy for each sample is shown in Table 5. 





Sample Total Enthalpy [eV] 





1024 Fe perfect crystal 


-4221.373163 





1024 Cr perfect crystal 


-3928.353757 





1023 Fe vacancy 


-4215.461833 





1023 Cr vacancy 


-3921.989661 





1025 self-interstitial <1 10> (Fe-Fe) in Iron 


-4221.974874 





1025 self-interstitial <1 10> (Fe-Fe) in Chromium 


-3929.176430 





1025 mixed interstitial <1 10> (Fe-Cr) in Iron 


-4222.064124 





1025 mixed interstitial <1 10> (Fe-Cr) in Chromium 


-3928.756038 





1025 self-interstitial <1 11> (Fe-Fe) in Iron 


-4221.546602 





1025 self-interstitial <1 11> (Fe-Fe) in Chromium 


-3929.176430 





1025 mixed interstitial <1 11> (Fe-Cr) in Iron 


-4221.896441 





1025 mixed interstitial <111>(Fe-Cr) in Chromium 


-3928.756244 





1025 self-interstitial <1 10> (Cr-Cr) in Chromium 


-3926.664855 








1025 self-interstitial <1 11> (Cr-Cr) in Chromium 





-3926.579530 





Table 5. Table of Final Values of Total Enthalpy for Pure Element Samples 
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The formation energies for the vacancy (Er ) in pure elements are then computed 
from the total enthalpy of the pure element sample C=) and the total enthalpy of the 
sample with the vacancy (Ee ), according to the following formula: 
pe 
1024 


The formation energies for self interstitials (Eg) were computed using a very 





i E* —1023-( ) (14) 


similar approach, but incorporating the fact that samples which contain the interstitial 
consist of 1025 instead of 1023 atoms, and therefore correlating the energies in the 


following way: 


int errstitia E pure 


interrstitial . 


where Fr ossaioms 18 the total enthalpy of the sample with 1025 atoms and a self interstitial 


and Ev =~ is the total enthalpy of the sample with the pure element in pure crystal 


1024atoms 
1024 atom arrangement. 
For the mixed interstitial samples, however, a different methodology must be 


applied, defining the formation energy of the interstitial as follows: 


1Cr int erstital 


EFI = Pisce: -_ | ee 


where 
reference = Freep cs from—-hof ( 1 6) 
and 
— : pureCr . pureFe 
interpolation —_ Nc, | ens a Ni | ee 
E from-hot » HOWEver, 1s a product of a polynomial fit to Ah[eV/atom] as seen in 


figure 3-4, and as we are working in the region below 6 % of Cr concentration, where 
this difference function is negative, we must only use the Ah values from a series of 
samples with an increasing Cr concentration, relaxed by running them through the Moldy 


code, as seen below: 
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xCr Ah[eV/atom] 

0 0 

0.00097 -9.65E-05 
0.006836 -0.000564056 
0.012695 -0.000492174 
0.015625 -0.00088443 1 
0.019531 -0.001196044 
0.030273 -0.000627474 
0.038086 -0.00066705 1 
0.047852 1.16076E-05 














Table 6. Values of Ah[eV/atom] as Obtained by Modeling 


When fitted with a 4" order polynomial using Excel, this produces the following heat of 


formation function in the sub-6% Cr concentration region: 


E-Eideal [mev/at] 


0.0002 





y = -200.11x* + 10.993x° + 1.6459x? - 0.0821x 
-0.0012 e 


Figure 21. Heat of Formation Fit Below 6% Cr concentration 


Therefore, when x is calculated as the value of x = 1/1025 and the corresponding y turns 


out to be -7.85209 X 10° meV/atom. E is the 1025" multiple of that value of y 


from—hof 
and is then entered into (16). 
The results of this characterization are displayed in Table 7. 
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Fe 1.72 3.52 3.95 3.43 3.61 - - 
Cr 2.56 2.61 2.61 3.23 3.23 5.52 5.61 


























Table 7. Formation Energies for the Pure Element Characterization 


C. EVALUATING THE VACANCY ENERGY OF FORMATION 


In order to analyze the effect of increased Cr concentration on the formation 
energy of a single vacancy, it was necessary to characterize this behavior with a series of 
samples of increasing Cr concentration. Also, since it is unknown what the influence of 
the distance between the solute Cr atoms and the vacancy is, it was necessary to analyze 
the entire ensemble of possible configurations. 

Therefore, 11 samples of 1024 atoms with Cr concentrations ranging from a 
singular Cr atom to 17 % were simulated using the molecular dynamics code Moldy [45], 
and set up for a constant temperature of 0.01 K, no pressure, no chemical potential, with 
free sample volume, and a run of 3000 steps with a time step of 10° seconds. These runs 
provided the reference values of total enthalpy for the pure crystal samples. Those values 


are shown in Table 8. 
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Cr Total Enthalpy 
Concentration | [eV] 
9.77E-04 -4221.45017 
0.01563 -4221.88415 
0.03027 -4221.0655 
0.04785 -4219.7705 
0.06055 -4217.80019 
0.07422 -421418631 
0.08984 -4210.24463 
0.09668 -4209.23258 
0.10645 -4206.00442 
0.11523 -4203.40938 
0.17188 -4184.4663 





Table 8. Total Enthalpy of the reference samples 


Thereafter, each of the samples was run through a sequence of two programs 
using a LINUX script. The first program went through the sample and removed one atom 
at a time, creating a vacancy. Only Fe atoms were removed in order to maintain the Cr 
composition. At the end of operation, approximately 1000 vacancy samples were ready, 
representing all possible configurations. Next, all of the newly created samples were 
simulated using the molecular dynamics code Moldy, and set up for a constant 
temperature of 0.01 K, no pressure, no chemical potential, with free sample volume, and 
a run of 3000 steps with a time step of 10° seconds. In the end, the LINUX code 
extracted the final values of total enthalpy from each sample output file. 

A histogram of energies could then be produced for each set of samples 


containing the same number of Cr atoms. Over this histogram, a Gaussian distribution 


function was fitted using ORIGIN™ with the canonical form of 





A wa) 


Y=Yot 
TW 
rth ea 
2 


ey) 





where yo, A, X- and w are fitting parameters. In our case, x. provides us with the mean 
value of the total enthalpy of the sample and w provides us with the width of the 2/3 
width of the Gaussian, which represents one standard deviation of results and gives us the 
width of the error bars. Each mean value for total enthalpy of a sample including a 
vacancy was then used in conjunction with the respective value of total enthalpy of the 


perfect crystal to provide the formation energy of the vacancy at the given Cr 


: 2 < : fi ‘ : : 1023 
concentration. Since this equation is first-order linear, the uncertainty inE as 


tot ? 
obtained from the Gaussian fit (w), transfers directly into Es. The vacancy formation 


energies obtained are shown in Table 9. 





























Cr Vacancy Energy of 
Concentration | Formation (eV) w/2 
9.77E-04 -1.713903803 | 0.00215 
0.01563 -1.715881926 | 0.0749 
0.03027 -1.71444999 | 0.01249 
0.04785 -1.722512512 | 0.01765 
0.06055 -1.761504721 | 0.05799 
0.07422 -421414418.6 | 0.06145 
0.08984 -1.785316452 | 0.06531 
0.09668 -1.799338475 | 0.06517 
0.10645 -1.838959668 | 0.06851 
0.11523 -1.812067625 | 0.06319 
0.17188 -1.863116579 | 0.05324 

















Table 9. Vacancy Formation Energy as a Function of Cr Concentration 


D. OBTAINING THE FORMATION ENERGIES OF INTERSTITIALS 


The procedure for obtaining the formation energies of self-interstitials is quite 
analogous to those described for the evolution of vacancy formation energy and for self- 
interstitial cases in the pure elements. After obtaining the reference values for the samples 
in the perfect crystal state, the samples are operated on by a program which places an 


interstitial atom added at each and every atom site containing a Fe atom. The program 
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displaces the original atom by 0.5 A away from its original location and adds its mirror 
image atom. The displacement takes place in x and y axes for interstitials <110> and all 
three axes for interstitials <111>. After creating approximately 1000 samples for each Cr 
concentration examined, the script then runs each of them in a Moldy code, with the 
inputs set at a constant temperature of 0.01 K, no pressure, no chemical potential, with 
free sample volume, and a run of 3000 steps with a time step of 10° seconds. The 
extracted values of final enthalpy are then sorted in an energy histogram and a Gaussian 
fit is applied, thus providing us with the mean energy value for each selected 
concentration as well as an indication of uncertainty based on the configurational 
characteristics of the sample. The values of final enthalpy for pure crystal samples, 
median values for those containing an interstitial, their individual uncertainties and the 


self-interstitial formation energy calculated via equation (15) are shown in Table 10. 





























Cr Formation Energy | w/2 Formation Energy | w/2 
Concentration | <110>eV <110> <111>eV <111> 
0.000976 3.519643423 | 0.004155 3.942243423 | 0.00729 
0.006829 3.516942656 | 0.006155 3.941462656 | 0.00883 
0.012683 3.512976597 | 0.019815 3.932626597 | 0.014195 
0.030244 3.506226279 | 0.031625 3.911206279 | 0.020885 
0.047805 3.504111963 0.04858 3.844751963 | 0.07725 
0.074146 3.494965576 | 0.106945 3.811195576 | 0.06379 
0.106341 3.498952173 | 0.154455 ; ; 
0.12 3.48620728 | 0.171505 ; : 
0.145366 3.4788 1083 0.19146 ; ; 
0.171707 3.47969771 | __ 0.19536 ; : 























Table 10. Energetics of Self-Interstitial Samples 


To obtain the evolution of mixed interstitial formation energies, we apply an 
almost exact procedure, except that when inserting an atom we are introducing a Cr, 
rather than a Fe. As this changes the global composition of the sample, we must correct 
for this using equations (16). In order to obtain Efrom por value it is necessary to make a 


dual zone polynomial fit. In this case zone A is below 6% Cr concentration and zone B 
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covers the area between 6-20% Cr concentrations. The relevant equations and fits are 


shown in Figures 22 and 23. 


0.002 5 y = 994.51x4 - 115.03x? + 6.0924x? - 0.1041x + 6E-06 


0.0015 + 
0.001 


0.0005 


delta h (eV) 








-0.0005 





-0.001 - 


Cr Fraction 


Figure 22. Polynomial Fit in the sub 6% Cr Concentration Eegion 
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Figure 23. Polynomial Fit to the 6-20 % Region of Cr Concentration 
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The resulting energetics of mixed interstitials are shown below in Table 11: 






































Cr Formation Energy Formation Energy 
Concentration | <110>eV w/2<110> | <111>eV w/2 <111> 
0.001951 3.428325907 0.007 3.596965907 0.004155 
0.007805 3.433025476 0.01334 3.591255476 0.021995 
0.013659 3.432935366 0.0165 3.603 105366 0.033825 
0.03122 3.467658628 0.06466 3.608058628 0.062045 
0.04878 3.509552851 0.0963 3.628022851 0.12394 
0.075122 3.550068635 0.105995 3.633258635 0.15496 
0.107317 3.546490688 0.10822 3.621460688 0.184035 
0.120976 3.537463852 0.109665 3.604573852 0.192875 
0.146341 3.530354706 0.126305 3.546284706 0.140005 
0.172683 3.484627439 0.124385 3.501687439 0.18874 

















Table 11. Mixed Interstitial Energetics in eV 
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VI. ANALYSIS AND CONCLUSION 


A. VACANCIES 


As we can see from Figure 24, our results for vacancy formation energy in Fe are 
lower than DFT calculations.. However, when we consider the fact that most of these 
calculations were either thermally unequilibrated (unrelaxed) or had a volume constraint, 
this is becomes a much more believable result. Both of these conditions led to additional 
“frustration” within the sample, and accordingly to a higher formational enthalpy. 
References for all values quoted are listed at the end of the chapter. 

On the other hand, when we compare the vacancy formation energy in Cr with 
other results (Figure 25), as well as those calculations which produce both results (Figure 
26), we see that it falls neatly in between ab initio and other classical results. From both 
of these comparisons we can conclude that our potential gives a reasonable prediction of 
vacancy formation energies for both pure elements, as well as for their mutual 
relationship. 

The more extensive analysis of vacancy formation energy behavior as a function 
of Cr concentration yields some fascinating results. Prior to this work it has been 
supposed that the vacancy formation energy would follow a direct interpolation between 
the pure element values. As we can see in Figure 27, the true formation energy deviates 
from the assumed linear interpolation in Figure 26 in the region below 6% Cr 
concentration, which is the same region in which the heat of formation curve shows the 
same deviant behavior from linear interpolation. In real terms, this finding means that it is 
more likely for vacancies to form in alloys with a Cr concentration smaller than 6%, then 
those above this value. The exact nature of this effect on the physical properties of the 
alloy is as of now undeterminable; however, it can be hypothesized that it will be 
observable. Additionally, the increasing spread of values as a function of Cr 


concentration shows a strong interactivity between the impurities and vacancies. 
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Figure 25. Vacancy Formation Energies in Chromium 
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Figure 26. Linear Interpolation Between EFV in Fe and Cr 
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Figure 27. Evolution of Vacancy Formation Energy as a Function of Cr 
45 


B. SELF-INTERSTITIALS 


As we can see from Figures 28 and 29 our results match both the relative order 
and the magnitude of target ab initio calculations. The critical part of this fitting was 
achieving the same relative order of values for mixed and self interstitials in iron, as both 
theory and experiment suggest the mixed interstitial to be the prime defect in Fe-Cr alloys 
[45]. 

Examining the plot of interstitial formation energy as a function of Cr 
concentration yields some exceptionally interesting information. While the formation 
energy of the interstitial in <110> remains relatively stable around 3.5 eV per defect, and 
actually drops lower as Cr concentration increases, the <111> curve shows some highly 
unstable behavior. An initial examination of energy histograms, produced by compiling 
and sorting the final output enthalpies of all samples initialized with a <111> interstitial, 
showed a behavior that deviated significantly from that which was expected. For 
example, Figure 30 shows the energy histogram for samples containing 10.6 % Cr atoms 
and an interstitial in <110>. The energies in the histogram follow a normal distribution 
which can be approximated and analyzed with a Gaussian curve fit as overlaid. For 
samples which contain just a single Cr atom the values are also normally distributed (see 
Figure 31); however, their variance is much smaller because most configurations interact 
equally with the interstitial. The only deviants are those where the Cr atom is in the near 
neighborhood of the interstitial site. 

However, when we take a look to the histogram of sample enthalpies with 0.006 
Chromium concentration, in Figure 32, we make the astonishing observation that there is 
not one, but two distinct distributions of energies. Upon closer examination, we notice 
that the larger distribution has a mean value which is equivalent to that of the distribution 
for samples with the <110> interstitial at the same Cr concentration. This inevitably leads 
us to conclude that a number of samples, even though initially seeded with a <111> 
interstitial, converted to an interstitial in the <110> direction during the process of energy 
minimization. Physically, this means that the <111> interstitial is not a preferred state for 
the defect and even if induced it will tend to convert into a <110> type interstitial. The 
fraction of samples which convert into a <110> type interstitial increases steadily as 
shown in Figure 33. This increase, combined with the fact that the differential between 
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self-interstitial formation energies in <110> and <111> decreases as seen in Figure 34, 
makes it impossible to discern the existence of a peak with final enthalpies for samples 
containing the <111> type interstitial after approximately 10 % because the variance has 


increased so much that the two distributions completely overlap. 


C, MIXED INTERSTITIALS 


When looking at the mixed interstitials and their formation energy evolution as a 
function of Cr concentration, we note several important differences when compared to 
the self-interstitials. 

In this case, it is the <110> interstitial that has a lower energy of formation, 
therefore being the preferable state. However, both mixed interstitial curves, as can be 
inferred from comparing Figure 35 with Figure 34, are above those for self-interstitials. 
However, they are both exceptionally stable and well behaved, so much so that they run 
almost parallel above 10 % in concentration and even have very little overlap when their 
variances are taken into account. 
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Figure 28. Formation Energies for Interstitials in Iron [45] 
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Figure 29. Formation Energies for Self-Interstitials in Chromium 
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Figure 30. Energy Histogram for Self-Interstitial Samples in the <110> direction at 


higher Cr Concentration 
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Figure 31. Histogram of Energies for Self-Interstitial Samples in the <110> direction at 


lower Cr Concentration 
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Figure 32. Histogram of Energies for Self-Interstitial Samples in <111> 
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Figure 33. Evolution Self-Interstitial Formation Energies 
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Figure 34. Partition Function for Stability of Self —Interstitial in <111> 
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Figure 35. Evolution of Mixed-Interstitial Formation Energies 
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Figure 36. Combined Interstitial Evolution Plot 
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D. CONCLUSIONS 


If we look at Figure 36, we can see that all interstitial types behave in a very 
similar manner higher the Cr concentration. It is not clear what physical significance this 
prediction could have, but it would have to be validated through either experimental 
methods or ab initio calculations, before further discussion can ensue. Nonetheless, the 
ability of EAM to offer such thorough evaluations should be noted. 

To conclude, it has been proven that this methodology [29] adequately 
approximates results obtained through ab initio calculations. It has also been 
demonstrated that by using this methodology, the studies of point defects can be of a 
larger scale and more exhaustive than is the case with the currently available ab initio 
models. 

The results produced in this research will go into studies of dynamic properties of 
point defects such as vacancy migration and defect clustering. In the meantime, though, a 
interim step is a study of precipitation properties of Fe-Cr alloys, which will have a high 
fidelity considering that the underlying assumptions, 1.e., the interatomic potentials, have 
already been validated as providing an accurate approximation of ab initio results, and by 
inference real world physical behavior. 

In order to study this behavior, a hybrid Monte Carlo-molecular dynamics code 
named MCCASK was developed by A. Caro and B. Sadigh in 2005 at Lawrence 
Livermore National Laboratory. The code performs sequences of Monte Carlo events and 
Molecular Dynamics time steps. In this way, the equilibrium concentrations in the alloy 
are obtained. 

The molecular dynamics part of the MCCASK code is based in MDCASK, a 
molecular dynamics code that was first developed to study radiation damage in metals. 
MDCASK solves the equation of motion of the atoms in the sample; energy and forces 
on each atom are calculated using an interatomic potential, and the equations of motion 
are integrated to obtain the next values of positions and velocities. The Monte Carlo part 
of the new code MCCASK corresponds to a parallel MC code in the transmutation 
ensemble (T,P,N,Dm), with N the total number of atoms and Dm the difference in 


chemical potential [46]. 
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Computational materials science is a young field, which has just recently been 
handed the tools of its trade--high-end computing hardware and software--and over the 
coming years it is inevitable that major advances in modeling will result in real life 
impact, not only in the field of nuclear energy, but in everyday applications. The future is 


bright. 
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VI. SUMMARY OF RESEARCH 


In summary, it can be said that the suitability of Fe-Cr alloys for use in cladding, 
wrappers and in-core structural components in Generation IV advanced nuclear systems 
necessitates a thorough understanding of their behavior under irradiation, to include 
formation energies of point defects such as vacancies and interstitials in the steel matrix. 
This work used a powerful methodology developed by A. Caro, et al. [29] to determine 
the energetics of vacancies and interstitials both in pure elements and as a function of Cr 


concentration. An overview of the methodology is shown in Figure 37. 
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Figure 37. Research Methodology [29] 
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Results of this research proved that the methodology is appropriate and can be 
utilized for large scale simulation, and that it approximates ab initio results satisfactorily. 
In addition to that, the evaluation of point defect formation energy as a function of Cr 
concentration, gave some interesting predictions concerning their static properties and 
stability. A finding relating the stability of the <111> Fe self-interstitial to local Cr 


concentration is summarized in Figure 38. 


5 - 
4.8 5 







4.6 5 
4.4 5 
4.2 + 


Stability Function of the <111> 
Interstitial 


f 


eRERRRERERE 











4 aa oe 
3.8 5 - 
3.6 - 
3.4 5 
52+ 
Se 
0 0.05 0.1 0.15 0.2 
xCr 


Figure 38. Development of the Stability Function for Fe Self-interstitial in <1 10> 


We can therefore say that the results obtained in this thesis can be used for studies 
of dynamic properties and time evolution of point defects. The stability predictions will 
need to be verified through ab initio or experimental data, but at the very least they can 
serve to demonstrate the potential benefits of using molecular dynamics and semi- 


empirical potentials to approximate first principles calculations. 
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